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ABSTRACT 

Relativistic astrophysical phenomena such as gamma-ray bursts (GRBs) and active galactic nuclei often re- 
quire long-lived strong magnetic field that cannot be achieved by shock compression alone. Here, we report 
on three-dimensional special-relativistic magnetohydrodynamic (MHD) simulations that we performed using 
a second-order Godunov-type conservative code, to explore the amplification and decay of macroscopic turbu- 
lence dynamo excited by the so-called Richtmyer-Meshkov instability (RMI; a Rayleigh-Taylor type instabil- 
ity). This instability is an inevitable outcome of interactions between shock and ambient density fluctuations. 
We find that the magnetic energy grows exponentially in a few eddy-turnover times, because of field-line 
stretching, and then, following the decay of kinetic turbulence, decays with a temporal power-law exponent 
of -0.7. The magnetic-energy fraction can reach ^ 0.1 but depends on the initial magnetic field strength, 
which can diversify the observed phenomena. We find that the magnetic energy grows by at least two orders 
of magnitude compared to the magnetic energy immediately behind the shock, provided the kinetic energy of 
turbulence injected by the RMI is larger than the magnetic energy. This minimum degree of the amplification 
does not depend on the amplitude of the initial density fluctuations, while the growth timescale and the maxi- 
mum magnetic energy depend on the degree of inhomogeneity in the density. The transition from Kolmogorov 
cascade to MHD critical balance cascade occurs at ^ 1/lOth the initial inhomogeneity scale, which limits the 
maximum synchrotron polarization to less than ^ 2%. We derive analytical formulas for these numerical re- 
sults and apply them to GRBs. New results include the avoidance of electron cooling with RMI turbulence, the 
turbulent photosphere model via RMI, and the shallow decay of the early afterglow from RMI. We also per- 
formed a simulation of freely decaying turbulence with relativistic velocity dispersion. We find that relativistic 
turbulence begins to decay much faster than one eddy-turnover time because of fast shock dissipation, which 
does not support the relativistic turbulence model by Narayan & Kumar 
Subject headings: magnetic fields — turbulence — instabilities — relativity — gamma rays 



1. INTRODUCTION 

It is well known that magnetic fields play a very important 
role in high-energy astrophysical phenomena such as parti- 
cle acceleration and synchrotron emission. Strong magnetic 
fields are often required around shock waves to explain the 
observed emissions. For example, at the shocks in super- 
nova remnants (SNRs), gamma-ray bursts (GRBs), and active 
galactic nuclei (AGNs), magnetic field strengths with orders 
of magnitude larger than their ambient fields are needed, and 
these cannot be achieved by only shock compression (e.g., 
Medvedev & Loeb 1999). Thus, some amplification mech- 
anisms have been proposed such as the Weibel instability 
(Weibel 1959), the cosmic -ray streaming instability (Lucek & 
Bell 2000; Bell & Lucek 2001), and the (small-scale) dynamo 
effect (Batchelor 1950) due to the turbulence induced by mag- 
netohydrodynamic (MHD) instabilities (Balsara et al. 2001; 
Giacalone & Jokipii 2007; Inoue et al. 2009, 2010; Zhang 
et al. 2009; Mizuno et al. 2010). Each mechanism has its 
own advantage. The Weibel instability can create magnetic 
fields even from unmagnetized media, and the cosmic -ray 
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Streaming instability can amplify upstream magnetic fields 
that are essential in first-order Fermi-acceleration processes 
in the shock wave. Because the dynamo effect works un- 
der the influence of macroscopic hydrodynamic instabilities, 
it can generate larger scale magnetic field fluctuations than 
those induced by the microscopic plasma instabilities, which 
may be essentially important for synchrotron radiation and 
scattering high-energy particles. Using three-dimensional rel- 
ativistic MHD simulations, Zhang et al. (2009) showed that 
the Kelvin-Helmholtz instability (KHI) induced by relativis- 
tic shear flows activates turbulence dynamos. 

In this study, we investigate the turbulent-dynamo effect 
with three-dimensional relativistic MHD simulations by tak- 
ing GRBs as an example. We study the effects of the 
Richtmyer-Meshkov instability (RMI), which are induced 
when the preshock density is inhomogeneous (see Brouillette 
2002; Nishihara et al. 2010, for reviews of RMI). The RMI 
is a Rayleigh-Taylor-type instability that deforms the shock 
front and leaves vorticity behind the shock waves. For the 
Rayleigh-Taylor instability, deformation of a discontinuity is 
triggered by gravitational acceleration, while the RMI is trig- 
gered by an impulsive acceleration caused by the passage of 
shock. 

Because GRBs are believed to be associated with rela- 
tivistic shocks of intermittent and inhomogeneous outflows 
(Meszaros 2006), we can reasonably expect the work of 
RMI. In the framework of Newtonian fluid dynamics, it was 
shown using MHD simulations in the context of SNRs that 
the preshock density inhomogeneity cause the magnetic field 
amplification (Giacalone & Jokipii 2007; Inoue et al. 2009, 
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2010). In a relativistic shock wave, the idea was examined 
analytically by Sironi & Goodman (2007), who considered 
how vorticity is generated because of the interaction between 
an ultra-relativistic shock and a density bump. Recently, us- 
ing two-dimensional special relativistic MHD simulations, 
Mizuno et al. (2010) showed that the magnetic field can be 
amplified in the postshock of an inhomogeneous medium. 
However, it is known that turbulence in two dimensions usu- 
ally leads to very intermittent structures of velocity and mag- 
netic fields compared to the three-dimensional case due to in- 
verse cascade of enstrophy (Biskamp 2008) that may induce 
"fake" magnetic field amplifications. Thus, the evolution of 
turbulence and resulting magnetic field amphfication should 
be verified by three-dimensional simulations, although the re- 
sult of Mizuno et al. (2010) is suggestive. 

In this paper, using three-dimensional special relativistic 
simulations, we examine the following items: (i) the gener- 
ation of turbulence at the relativistic shock front through the 
RMI and the resulting magnetic field amplification caused by 
turbulence and (ii) the decay of relativistic MHD turbulence 
based on the results of (i). Because it is computationally very 
expensive to follow the long-term evolution of turbulence in 
the former simulations, we cannot obtain the saturation level 
or decay rate of the turbulent magnetic field. The supplemen- 
tary simulations of (ii) allow us to study the decaying phase 
of turbulence as a consequence of its long-term evolution. 

This paper is organized as follows. In §2 we provide nu- 
merical techniques and setups. The results of the numerical 
simulations are presented in §3 (generation of turbulence and 
growth of magnetic field) and §4 (long-term evolution and 
turbulence decay). In §5, we summarize our findings from 
the simulations. Finally, in §6, we discuss the implications of 
the GRB phenomenologies, including the internal shock syn- 
chrotron model, the photosphere model, the afterglow, and 
polarizations of these phenomena. 

2. NUMERICAL SETUP 

We solve the ideal special-relativistic MHD equations. Be- 
cause we treat propagation of relativistic shock waves and 
turbulent magnetic fields, a high-resolution shock-capturing 
scheme as well as a divergence-free induction-equation solver 
are preferred. We have developed such a code by using the 
five-wave Harten-Lax-van Leer Riemann solver (Miyoshi & 
Kusano 2005) developed by Mignone et al. (2009) and the 
constrained transport scheme (Evans & Hawley 1988) de- 
signed by Stone & Gardiner (2009). The five- wave HLL 
solver employed in this paper is also called the HLLD solver. 
According to Beckwith & Stone (201 1), HLLD solver is more 
suited to treat magnetic field amplification by turbulence than 
other HLL solvers. We impose the TM equation of state 
(EOS) proposed by Mignone et al. (2005), which correctly 
describes the EOS for nonrelativistic and relativistic temper- 
ature limits and can describe the transition regime that differs 
from the exact EOS by less than 4%. We solve the equations 
in the conservative fashion so that the mass, momentum, and 
energy are conserved to within the round-off error. Because 
we solve the scale-free, ideal (relativistic) MHD equations, 
we use a dimensionless time normalized by the light-crossing- 
time of the numerical domain t = tc/L,, where L- is the length 
of the numerical box along the z-axis. Thus, for example, if 
we choose the scale = 10^ cm, f = 1 corresponds to 0.33 
msec. For intuitive presentation, we treat variables such as 
density, velocity, pressure, and magnetic field dimensionally. 
Note that the choice of does not change the value of the 



TABLE 1 
Model list 



IxUll IlalllC 


Density dispersion 




Al 


An/no = 0.9 


7.50 G 


A2 


0.6 


7.50 G 


A3 


0.3 


7.50 G 


A4 


0.1 


7.50 G 


Run name 


Velocity dispersion* 


Magnetic field strength'^ 


Bl 


Av/cs = 1.0 


2170 G 


B2 


1.0 


217 G 


B3 


1.0 


21.7 G 


B2-S 


0.3 


217 G 


B2-r 


1.7 ((r) = 3.0) 


217 G 



^Initial magnetic field strength in the fluid rest frame of the colliding flows. 

''Cs = 0.48 c is the speed of sound in the initial medium. 

'^Initial magnetic field strength in the simulation frame, i.e., the post- 
shock rest-frame when upstream is uniform. Thus, the initial magnetic field 
strengths shown correspond to the values obtained after the compression by 
internal shock. 

above variables-it only affects the timescale. 

2.1. Initial Conditions of the First Set of Simulations 

In the first set of simulations, we prepare a cubic domain 
whose aspect ratio is : Lj, : Lj. = 2 : 1 : 1, which is divided 
by iVx X iVy X A^z = 5 12 X 256 x 256 finite uniform volume ele- 
ments. We initially generate an inhomogeneous medium by 
adding isotropic sinusoidal density fluctuations in the fluid 
rest frame: 

n = no+ ^ P{k)^^^ ?.m{kxX+kyy+k^z+(t>k), (1) 

where (j)k is a random phase that depends on the wave num- 
ber. The mean number density of the medium is set to 
(n) = no = 10^" protons cm~^, which is compatible with GRB 
shells that collides at r ^ lO'"* cm from the central engine. 
The fluctuations are characterized by their power spectrum 
Pik). We fix the spectral shape as a flat spectrum k^P{k) oc k^ 
cutoff at the scale 



so that the fluctuations are only in large scales. Note that the 
suimnations are performed for nonzero k = {k]. + k^ + k\)^/'^. 
We examine four situations that have density dispersions 
An/(n) =0.1, 0.3, 0.6, and 0.9, where An = ^{n^)-{n)^. 

To induce shock waves, we set the converging velocity field 
to = 0.8 1 8 c at X < Lj./2 and v.j. = -0.8 1 8 c at x > Lj./2, where 
c is the speed of light, which indicates that the simulations 
are performed in a postshock rest frame when the upstream 
is uniform. The Lorentz factor of the flows in the simulation 
frame is F = 1 .74 (the relative Lorentz factor of the flows is 
r = 5.05), which corresponds, for example, to the collision 
of the GRB shells with Liab = 1000 and 100 in the laboratory 
frame. 

The pre- shock temperature of the GRB shell in the fire- 
ball model is usually too cold for our conservative numeri- 
cal scheme to perform stable simulations (see eq. [9]). Thus, 
we increases the temperature to 9.38 MeV for the purpose of 
numerical stability. By solving shock-jump conditions, we 
have confirmed that the discrepancy of the average postshock 
conditions (density, pressure, and magnetic field) between the 
case of this warm preshock temperature and the cold-limit 
case (r = K) is only 3% because our initial temperature is 
still much lesser than the postshock temperature of ~ 1 GeV. 
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Thus, we stress that our initial condition is sufficiently cold to 
generate realistic GRB internal shock propagation. 

The initial magnetic field strength in the fluid rest frame is 
fixed at B = 7.50 G (B = 13.0 G in the simulation frame) ori- 
ented in the +y direction. As discussed in §6, if we assume the 
frozen-in transport of magnetic field from the central engine, 
this initial strength corresponds to a central engine magnetic 
field of ^ 4 X lO'" G (see, eq. [8]). Periodic boundary condi- 
tions are imposed for the y and z boundaries. At the x bound- 
ary, we continue to impose the converging relativistic flows 
in which density fluctuations same as the initial fluctuations 
are periodically input. The initial model parameters and the 
names of the runs are summarized in Table 1 (see also §6 for 
the corresponding GRB initial conditions). 

2.2. Initial Conditions of the Second Set of Simulations 

As shown in the next section, the set of simulations dis- 
cussed in the preceding section can only treat the generation 
of turbulence and the growing phase of the magnetic field, 
because the induced shock waves rapidly exit the computa- 
tional domain. To study the decay of turbulence and magnetic 
field, we need to perform long-term simulations. To do so, 
we perform the following set of simulations: We prepare a 
cubic domain with : Ly : = \ : \ : \ having resolution 
X Ny X Nz = 512^. We fill the domain with an uniform 
medium whose density and temperature are « = 10'' cm"-' and 
k^T = 0.235 GeV, respectively, which are compatible with the 
averaged postshock values obtained from the first set of sim- 
ulations. Turbulence is initially given by summing the sinu- 
soidal velocity fluctuations with zero mean velocity: 



rv,= 



E 



P{ky^^ sin{k.tX+kyy+k^z+(f>k.d, 



(3) 



where the subscript / covers x, y, and z, and (f)^., is a random 
phase that depends on the wave number and the components. 
Note that the summation is overall nonzero A: = (fc^+A;^+A;?)'''^, 
and we generate fluctuating fields of F v, to avoid creating su- 
perluminal regions. Their power spectrum is the flat with cut- 
off scale 

In the generated initial velocity field thus generated, the ratio 
of the total power of the rotational velocity component ( V • v*= 
0) to the compressive component ( V x v = 0) is approximately 
2:1, which indicates that the initial turbulence is dominated 
by incompressible flows. The periodic boundary conditions 
are imposed and we perform several runs with different initial 
velocity dispersions (Av = \/Jy^) and initial magnetic field 
strengths, whose values are summarized in Table 1 (see also 
§6 for the corresponding GRB initial conditions). 

3. SHOCK PROPAGATION IN INHOMOGENEOUS MEDIA 

3.1. Generation of Turbulence 

The initial converging flows used in the series of run A 
induce two oppositely propagating shock waves. When the 
shock wave passes a density bump or dent, a RMI is in- 
duced and leaves vorticity behind the shock waves. For the 
Rayleigh-Taylor instability, deformation of a discontinuity 
and generation of vorticity are triggered by gravitational ac- 
celeration, while the RMI is triggered by an impulsive accel- 
eration caused by the passage of the shock wave (Brouillette 
2002; Nishihara et al. 2010). 




^1 



0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 




Fig. 1. — (a): Two-dimensional number density slice of the result of run 
Al. (b): Two-dimensional magnetic field strength slice of the result of run 
Al. (c): Two-dimensional number density slice of the result of run A3, (d): 
Two-dimensional magnetic field strength slice of the result of run A3. All 
panels show the structures al z = 0.0 and t = 2.5. Physical quantities are 
measured in the simulation frame, which is equivalent to the postshock rest 
frame when the upstream is uniform. 

Fig. 1 (a) and (c) show two-dimensional number-density 
slices of the results of runs Al and A3, in which one can see 
the deformed shock fronts and turbulent density fields due to 
the RMI. In the following, we plot the values measured in the 
simulation frame in figures, which is equivalent to the post- 
shock rest frame when the upstream is uniform. Since lower 
(higher) upstream density leads to faster (slower) shock prop- 
agation, the medium with larger density dispersion induces 
larger shock deformation, and thus leaves stronger vorticity or 
turbulence. To clarify this, we plot the dispersion of the sonic 
Mach number in the y-z plane at a given x in Fig. 2. The red, 
blue, green, and magenta lines show the results of runs Al, 
A2, A3, and A4, respectively. Series of the same color lines 
represent different temporal snapshots at f = 0.5, 1 .5, and 2.5. 
We also plot the reference lines of M = 0.1, 0.3, 0.6, and 0.9 
as dashed lines. It is obvious that larger preshock density 
dispersion introduces stronger postshock velocity dispersion, 
and that the velocity dispersion is roughly proportional to the 
preshock density dispersion. As the postshock Mach number 
dispersion approaches unity, in particular run Al, the Mach 
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Fig. 2. — Dispersion of the sonic Mach number in the y-z plane at a given 
X. The red, blue, green, and magenta lines represent the results of runs Al, 
A2, A3, and A4, respectively. Series of the same-color lines show different 
temporal snapshots at ? = 0.5, 1.5, and 2.5. Also plotted as dashed lines are 
the reference lines for M = 0.1, 0.3, 0.6, and 0.9. 

number dispersion becomes somewhat smaller than that ex- 
pected from the linear dependence on the density dispersion. 
This would be due to the formation of "eddy shocklets" in the 
turbulent postshock (Kida & Orszag 1990) that tend to keep 
the velocity dispersion subsonic. 

Sironi & Goodman (2007) developed an analytic formula 
for the generation of vorticity because of an interaction be- 
tween a relativistic shock and a density bump using the geo- 
metric shock-dynamics approximation (see, also Goodman & 
Macfadyen 2008). They found that the vortical energy den- 
sity generated by the interaction is proportional to the square 
of the maximum density of the bump (in other words, the 
velocity dispersion of vorticity is proportional to the den- 
sity contrast), when the density contrast is small. Because 
they assumed ultra relativistic shock wave to study the af- 
terglow phase of GRBs, it may be inappropriate to compare 
their result and that of the mildly relativistic case examined 
here. However, note that our result (linear dependence of the 
velocity dispersion on the preshock density dispersion when 
An/no <C 1) is consistent with their finding. In addition to the 
density bump-shock interaction, also note that the interaction 
with density dents also generates vorticity because the RMI is 
known to be effective even in that case. 

Fig. 2 shows that even when f = 2.5 the level of the post 
shock velocity dispersion at around x/L, = 1.0 is almost the 
same as for f = 0.5, which indicates that the turbulence has not 
yet begun to decay. In the present simulations, the postshock 
velocity dispersion is nonrelativistic, and in Newtonian fluid 
dynamics, it is widely known that the turbulence without a 
continuous driving source begins to decay within a few eddy- 
turnover times. Because the average postshock sound speed is 
Cs ~ 0.5 c and the scale of vortices that is given essentially by 
the scale of the density fluctuations is ^ L., the eddy-turnover 
time can be estimated to be a few light-crossing times of L.. 
Thus, to confirm the decay of turbulence, we need to continue 
the simulation for several light-crossing times, which requires 
a larger numerical domain. The decaying phase of turbulence 
will be shown in §4 as the results of the second set of simula- 
tions. 

3.2. Magnetic Field Amplification 

In addition to the shock compression, the turbulence in- 
duced by the RMI amplifies the magnetic field by stretching 
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time : t 

Fig. 3. — Top: Evolution of the maximum (|B|max: lines) and average 
((|B|): points) magnetic field strengths. Bottom: Evolution of average mag- 
netic energy density at the x/L- = 1.0 plane ((r^B^/87r),_i). Red, blue, 
green, and magenta lines show the results of runs Al, A2, A3, and A4, re- 
spectively. Dashed lines show the model evolution eq. (5) with / = 1.8. 

the field lines. The slices of magnetic field strength distri- 
bution of runs Al and A3 are shown in Fig. 1 (b) and (d), 
respectively. We can confirm the magnetic field amplifica- 
tion far beyond the value achieved by the shock compression 
alone (if the preshock medium is uniform, the magnetic field 
strength is B 70 G). 

Similar amplification of magnetic fields in the dynamics of 
SNRs formed in various ambient ISM has been reported by 
many authors (see Balsara et al. 2001 for turbulent ISM; 
Giacalone & Jokipii 2007 for ISM whose density fluctua- 
tion power spectrum is a power-law, which is recently ap- 
plied to relativistic shocks by Mizuno et al. 2010; Inoue et 
al. 2009, 2010 for cloudy ISM). The initial density fluctu- 
ations for our simulations are similar to that of Giacalone 
& Jokipii (2007) and Mizuno et al. (2010). However, they 
examined shock propagation in a two-dimensional geometry, 
and it is widely known that a two-dimensional turbulence can 
be qualitatively different from a three-dimensional turbulence 
because of the inverse cascade of enstrophy, which creates 
long-lived, large-scale eddies. At later stages, the large-scale 
eddies decay very intermittently and leads to a "fake" mag- 
netic field amplification (see, e.g., Biskamp 2008). Thus, the 
amplification reported in the previous two-dimensional sim- 
ulations seems suspicious (Beresnyak et al. 2009). Never- 
theless, such a dimensional difference stands out in the later- 
stage of driven turbulence, whereas the amplifications caused 
by density fluctuation-shock interaction is a short-term phe- 
nomenon. The presence of the magnetic field amplifications 
in our three-dimensional simulations suggests that the am- 
plifications reported in the previous two-dimensional simula- 
tions are not fake, although the results of the two-dimensional 
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simulations may be affected by the intemittency in some de- 
gree. Recently, Inoue et al. (2009, 2010) also showed that 
similar magnetic field amplification caused by shock-cloud 
interactions does not depend on spatial dimensions. 

We plot the evolution of the maximum and average mag- 
netic field strengths (|B|max and in the top panel of 
Fig. 3 and that of average magnetic energy density at the 
jc/Lj. = 1.0 plane ((T^ /%tt)x=i) in the bottom panel, where 
r is the Lorentz factor of each fluid element. Note that, since 
the post shock turbulence is subsonic and thus nonrelativis- 
tic, r in the above definition of the magnetic energy den- 
sity is approximately unity. It is clear that the amplifica- 
tion still continues, since the turbulence that is the source 
of the amplification has not yet decayed. Since the mag- 
netic field is passive with respect to the turbulent velocity 
field, we can reasonably expect from the induction equation 
[dB/dt = V X (v X B) ~ AvB/Zturb] that the magnetic energy 
density evolves as 



(a) t=0.60 



t= 1.00 



{es) 



7Av 
i exp I 1 



tmb 



(5) 



where (eB)ini = (r^B^)ini is the initial magnetic energy den- 
sity (immediately behind the shock), /tmb is the scale of the 
turbulent flow, Av is the velocity dispersion of the turbulence, 
and / is a factor on the order of unity. In the present simula- 
tions, the turbulent flows with the scale /mrb = L^/^ would give 
the fastest growing mode, because the scale of the flows is es- 
sentially determined by the scale of the initial density fluctua- 
tions. From Fig. 2, the velocity dispersion can be evaluated as 
Av/cs - 0.9, 0.6, 0.3, and 0.1 for the runs Al, A2, A3, and 
A4, respectively. In the bottom panel of Fig. 3, we plot the 
evolution of this simple model with / = 1.8 as dashed lines. 
We see that the model equation reproduces the growth rate es- 
pecially for f ^ 1 . In the next section, we show that the same 
model with / = 1.8 also reproduces the initial growth phase of 
magnetic energy, even for the second set of simulations. 

4. DECAY OF TURBULENCE AND MAGNETIC FIELD 

In the previous section, we have seen that the density 
fluctuation-shock interactions generate turbulence whose ve- 
locity dispersion can be as large as post shock sound speed 
depending on the amplitude of the density fluctuations. In this 
section, we show the results of the second set of simulation. 
The initial conditions of the run Bl, B2, B3, and B2-s, which 
have transonic and subsonic velocity dispersions, simulate 
the post shock turbulence induced by the density fluctuation- 
shock interactions. We can study the long-term evolution of 
turbulence and magnetic field from these simulations. 

The run B2-r has a relativistic initial velocity dispersion 
((r) = 3.0) that may not be created by the density fluctuation- 
shock interactions. However, since there were several dis- 
cussions on the possibility of relativistic turbulence in GRB 
(Lyutikov 2006; Narayan & Kumar 2009; Kumar & Narayan 
2009; Lazar et al. 2009), it would be meaningful to study it. 

4. 1 . Transonic and Subsonic Turbulence 

4.1.1. Evolution Process 

As in the first set of simulations, the magnetic field in the 
turbulence is amplified by the field-line stretching. In Fig. 4, 
we plot the z = plane slices of the magnetic field strength 
for run Bl. Panels (a), (b), (c), and (d) represent tempo- 
ral snapshots at F = 0.6, 1.0,4.0, and 8.0, respectively. The 
evolution of the maximum and the average magnetic field 




Fig. 4. — Slices of magnetic field strength of run Bl in the z = 0.0 plane. 
Panels (a), (b), (c), and (d) represent temporal snapshots at F = 0.6, 1.0, 4.0, 
and 8.0, respectively. 
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Fig. 5. — Top panel: evolution of the maximum and the average magnetic 
field strengths (|B|„,ax and obtained from runs Bl(red), B2 (blue), B3 

(green), and B2-s (magenta). Solid lines and points respectively indicate the 
maximum and average field strengths, respectively. Bottom panel: evolution 
of magnetic (solid), kinetic (dashed) and internal (dotted) energy densities. 
Line colors identify the models as per the top panel. Dashed black lines show 
the model of the evolution of the magnetic energy density given by eq. (5) 
with / = 1.8. The thin black line is a reference Hne proportional to f"" ' fit to 
the decay-phase evolution of the magnetic energy, and the dot-dashed black 
line is a reference line proportional to fit to the decay-phase evolution of 
the kinetic energy. 
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Fig. 6. — Evolution of eg obtained from runs Bl{red), B2 (blue), B3 
(green), and B2-s (magenta). 

Strengths (|B|max and obtained from runs Bl(red), B2 

(blue), B3 (green), and B2-s (magenta) is plotted in the top 
panel of Fig. 5. Solid lines and points indicate the maximum 
and average field strengths, respectively. In the bottom panel, 
we also plot the magnetic (solid), kinetic (dashed) and inter- 
nal (dotted) energy densities, where we have defined them 
(r^B^/STr), (pc2 r(r- 1)>, and ({F^ 7/(7- 1)- 1} respec- 
tively. Note that their summation is conserved to within the 
round-off error owing to the conservative numerical scheme. 
For runs Bl, B2, and B3, the kinetic energy begins to decay 
at f ^ 2 (i.e., a few eddy-turnover times), and the magnetic 
energy also begins to decay when it becomes comparable to 
the decaying kinetic energy. The initial velocity dispersion of 
run B2-S is approximately one-third of runs Bl, B2, and B3, 
which results in a larger eddy-turnover time and thus kinetic 
energy decay postponed by a factor of approximately three. 
Because magnetic energy does not dominate energy budget, 
results with different initial magnetic field strengths (i.e., runs 
Bl, B2, and B3) show almost the same evolution for the ki- 
netic and internal energies. 

Until the turbulence begins to decay, eq. (5) with / = 1.8 
(where we have substituted the initial velocity dispersion as 
Av and /turb = L^J 4- from the initial condition), which is plotted 
in Fig. 5 with a thin dashed line, again gives a suitable model. 
Note that in the case of a turbulence dynamo with a continu- 
ous large-scale driving source, the linear growth stage of the 
magnetic energy is followed by the initial exponential growth 
stage that continues until the magnetic energy becomes com- 
parable to the kinetic energy (Schekochihin & Cowley 2007; 
Cho et al. 2009). In our simulations, the turbulence is injected 
only initially and decays eventually, which halts the magnetic 
field amplification and leads to damping. In the decay phase, 
the evolution of the magnetic energy can be fit by a power 
law with an index ^ -0.7 and the kinetic energy also shows 
power-law decay with an index ^-1.3. The ratio of magnetic 
energy to internal energy, which is often denoted as eg is an 
important parameter in GRB-emission models. We plot the eg 
for runs Bl, B2, B3, and B2-s in Fig. 6. Because the initial 
internal energy is comparable to or larger than the initial ki- 
netic energy, the internal energy increases only slightly, even 
after the decay of turbulence. Thus, the evolution of the 
mainly determined by the evolution of the magnetic energy, 
and the simple growth model obtained by dividing eq. (5) 
by the initial internal energy well fit the early-phase evolu- 
tion (see, dashed lines in Fig. 6). The later power-law decay 
with the exponent -0.7 (dotted line in Fig. 6) also fits the later 
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Fig. 7. — Power spectra of velocity Pr{k) (top) and magnetic field Pgik) 
(bottom) at f = 0.2, 0.6, 1.0, 4.0, and 8.0 from run B2. Here, the power spec- 
tra are defined as / P,.(k)dk = Jv-d^x and / PB(k)dk = f B^d^x. Dashed 
Hne shows the Kolmogorov spectrum (Pv(k) oc k^^/^). 

evolution of eg- 

Can we predict the maximum eg using parameters such as 
Av and e^ini? In the present simulations, we can describe 
the eddy-turnover time for the smallest initial eddy as feddy — 
L,/4 Av (~ 0.5 Ljc for runs Bl, B2, and B3, and ~ 1.7L,/c 
for run B2-s, where we have used ~ 0.5 c). Thus, as seen 
from the bottom panel of Fig. 5, the turbulence can maintain 
its the initial strength until t ~ 3 feddy Substituting the above 
timescale into eq. (5), we obtain 



(eB(fdec)> 



:exp(3/)^102 



(6) 



Note that this degree of amplification describes the ratio of 
the magnetic energy at f = (immediately behind the shock 
wave) and at the time when the turbulence begins to decay 
(fdec)- This degree of amplification is independent of the initial 
velocity dispersion that would explain the comparable maxi- 
mum es of run B2 and B2-s. However, this degree of amplifi- 
cation cannot be used for the maximum eg prediction because 
the eb can grow, even after the fdec, until the magnetic en- 
ergy becomes comparable to the kinetic energy, although after 
the fdec the growth rate becomes slower than the exponential 
growth given by eq. (5). Indeed, for run B3, the maximum 
eg is approximately three orders of magnitude larger than the 
initial value. An accurate description of the growth of for 
t > fdec may not be obtained in a straightforward way; nev- 
ertheless, we can predict the upper limit of the maximum eg 
as follows: Because the growth of eg stops when the mag- 
netic energy becomes comparable to the kinetic energy, the 
upper-limit of should be smaller than the ratio of the ki- 
netic energy to the internal energy eK. For t > fdec, the kinetic 
energy ck evolves as oc f"'^, which leads to the inequality. 



< 



eK~ (M>ini(l+f/fdec)' 



-1.3 



(7) 



where (M)i„i is the Mach number of the initial turbulence (see 
Fig. 5). 



3D MHD TURBULENCE BEHIND RELATIVISTIC SHOCK WAVES 



7 




10' IC 
magnetic field strength [ Gauss ] 



10= 



Fig. 8. — Probability distribution histograms of the run Bl at ? = 0.4 (red), 
0.8 (blue), 1.2 (green), 1.6 (magenta), and t > 2.0 with the interval AF = 0.4 
(thin black lines). We also plot a distribution oc exp(— B/3 X 10^) (light- 
blue dashed line) as a reference. 
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Fig. 9. — Structure function defined by (|fo(r)- fo(r + /)p)f calculated by 
using the result of the run Bl at F = 4, where h = B/\B\. The vertical axis 
/|| indicates the distance parallel to the local magnetic field from r, and the 
horizontal axis /x indicates the perpendicular distance. 

4.1.2. Structure of Magnetic Field 

The spectra of velocity and magnetic field also evolve 
with time. In Fig. 7, we plot the power spectra of the ve- 
locity Pv{k) {top) and the magnetic field Pgik) {bottom) for 
t = 0.2,0.6, 1.0,4.0, and 8.0 of run B2, where the power 
spectra are defined as J Pv{k)dk = J v^dx^ and J PB{k)dk = 
J B^dx^ . The amplitude of the spectra evolves along with the 
kinetic and magnetic energy densities, whereas their shapes 
are roughly maintained after f ^ 1 (after a few eddy-turnover 
times). The shapes of spectra at f > 1 resemble the steady 
spectra of super-Alfvenic turbulence with a large-scale con- 
tinuous driving source (e.g., Cho & Vishniac 2000; Cho et al. 
2009; Zhang et al. 2009). For hydrodynamics, the spectrum 
of the rotational velocity component of the developed, decay- 
ing turbulence decreases its power while maintaining the Kol- 
mogorov spectrum (Kida & Orszag 1992). Recall that, in the 
series of run B, the initial power of the velocity field is dom- 
inated by the rotational component. Thus, it is not surpris- 
ing for the present spectra to evolve while maintaining their 
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Fig. 10. — Power spectra of velocity Pv{k) (top) and magnetic field Psik) 
{bottom) of run Al at ( = 2.5. Dashed Une shows the Kolmogorov spectrum. 

shapes. The spectra of other models of run B also evolve in a 
qualitatively similar way. 

In large scales {k/li: < 30/L,), the velocity power spec- 
trum for f > 1 exhibits a power law whose exponent is con- 
sistent with the Kolmogorov index of -5/3. However, the 
velocity power spectrum becomes steeper roughly at the scale 
~ L-/3Q, even though it is much larger than the scale of nu- 
merical resolution ~ L /500. This would indicate the transi- 
tion of the kinetic energy transfer mechanism from the hydro- 
dynamic Kolmogorov cascade to the magnetohydrodynamic 
critical balance cascade (Goldreich & Sridhar 1995). The 
power spectrum of the magnetic field also changes its slope 
at the transition scale. For large scales {k/ln < 30 /L,) Pb for 
f > 1 is roughly flat as commonly seen in the results of turbu- 
lence dynamo simulations (Cho & Vishniac 2000; Cho et al. 
2009; Zhang et al. 2009; Brandenburg & Subramanian 2005, 
and references therein), and at small scales it becomes steeper 
The steepening of Pg below the transition scale indicates that 
the magnetic energy is comparable to the kinetic energy at that 
scale. Appearance of the transition scale is quite reasonable, 
because the velocity dispersion of turbulent eddies decreases 
as the eddies cascade to smaller scales and eventually van- 
ishes, whereas magnetic field strength is not (i.e., there always 
exist a transition scale at which the velocity dispersion of ed- 
dies is comparable to the Alfven velocity constructed using 
the strength of the mean magnetic field). Thus, even if ve- 
locity dispersion is super-Alfvenic at large scales, it becomes 
sub-Alfvenic below the transition scale. 

The structure in the magnetic field strength in Fig. 4 sug- 
gests that the field strength distribution has large dispersion. 
In Fig. 8, we plot the probability distribution histograms of 
the magnetic field strength calculated using the result of run 
Bl at r = 0.4 (red), 0.8 (blue), 1.2 (green), 1.6 (magenta), and 
t > 2.0 with the interval Af = 0.4 (thin black lines). For f > 1, 
the distribution function can be fitted by the following func- 
tion: f{B) oc exp{-B/a{t)), where a{t) is the function of 
time. By taking the second moment of the normalized dis- 



8 



T. INOUE ET AL. 



tribution f(B) to be proportional to the decaying-phase mag- 
netic energy density (oc f~°'^), we obtain a{t) oc f"°-^^, which 
indicates that both the width and the peak of distribution de- 
creases with time for f > 1 as expected. We confirmed that the 
magnetic field strength distributions of the results of the other 
simulations (series of runs A and B) exhibit distributions sim- 
ilar to f(B) given above. 

Because the magnetic field is passive with respect to large- 
scale turbulent flows, the angular distribution of the magnetic 
field is nearly isotropic, whereas it shows a only slight de- 
pendence on the initial direction. However, the isotropic dis- 
tribution does not indicate the absence of local correlations 
of the magnetic field direction. In Fig. 9, we show the struc- 
ture function defined by {\b{r)- b(r+ /) p) ? calculated by using 
the result of run Bl for f = 4, where b = B/\B\. The vertical 
axis /|| gives the distance from r parallel to the local magnetic 
field, and the horizontal axis l± gives the distance perpendic- 
ular to the local magnetic field. If magnetic field orientation 
has spatial correlations, the structure function is null. Fig. 9 
shows that the magnetic field orientations have correlations 
on the small scale because the magnetic field is not passive 
below the transition scale. The anisotropic structure along the 
local magnetic field direction is qualitatively consistent with 
the theory of anisotropic cascade of Alfven waves (Goldreich 
& Sridhar 1995) and the result of turbulence dynamo simula- 
tions by Cho & Vishniac (2000). 

Finally, we note here that the spectra from run Al at f = 
2.5 that are plotted in Fig. 10 are similar to those of run B2 
at F > 1 in Fig. 7. This similarity suggests that the second 
set of simulations (i.e., the series of run B) is appropriate for 
studying the long-term evolution of the turbulence driven by 
the RMI. 

4.2. Relativistic Turbulence 



(a) 



For run B2-r, which has a relativistic initial velocity dis- 
persion ((r) = 3.0), the magnetic field also grows with time. 
We plot the z = plane slices of the magnetic field strength 
and number density for run B2-r at f = 0.10 in Fig. 11. In 
Fig. 12, we show the evolution of the maximum and average 
magnetic field strengths {top), the kinetic, magnetic, and in- 
ternal energy densities {middle), and {bottom). In run B2-r, 
the simulation is terminates at f = 0.16 (slightly after the ki- 
netic energy begins to decrease). The reason for this is as fol- 
lows: Initially, we set the relativistic turbulent field in which 
the regions with oppositely oriented relativistic flows that tend 
to create vacua there always exist. Because the employed nu- 
merical scheme is a conservative grid-based method, it is very 
difficult to treat a thin medium. The continuation of the sim- 
ulation might be possible, if we artificially put the mass and 
internal energy into the thin regions. However, this strategy 
may significantly changes the dynamics because the volume- 
filling factor of such thin regions can be large. Thus, we do 
not choose the continuation. 

Although the timescale we have followed is short, the ki- 
netic energy has already started to decay, whereas the mag- 
netic energy is still in the growth phase. The model equation 
(5) with / = 1.8, Av = c, and Zturb = L^/A, which is plotted 
as a dashed line in the middle panel of Fig. 12, again shows 
good agreement until the onset of the kinetic energy decay. In 
the present case, the decay of the kinetic energy begins ap- 
proximately ten times faster than that for the transonic cases, 
which we attribute to the shock dissipation that immediately 
converts kinetic energy into internal energy. Indeed, a number 
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Fig . 11. — Two-dimensional slices at z = of magnetic field strength struc- 
ture (top) and number density structure (bottom) from run B2-r at F = 0.10. 
Physical quantities are measured in the simulation frame. 

of discontinuous density structures shown in Fig. 1 1 indicate 
the formation of multiple shock waves. The beginning of the 
kinetic energy decay much faster than the eddy-turnover time 
is consistent with the simulations of decaying supersonic tur- 
bulence (e.g., MacLow et al. 1998), although these are non- 
relativistic, isothermal simulations. Note that, when the ki- 
netic energy injected into the turbulence is fixed, the timescale 
of the decay of the relativistic turbulence would substantially 
shorter than the nonrelativistic one, since the relativistic effect 
reduces the necessary mass for deceleration by F^. 

The evolution of eg is completely different from the tran- 
sonic and subsonic cases (i.e., €b decreases with time) be- 
cause the shock compression increases the internal energy 
more rapidly than the magnetic energy. As for transonic and 
subsonic turbulence, the magnetic energy would be able to 
grow until it becomes comparable to the kinetic energy. Thus, 
es can be turn to increase once the velocity dispersion of tur- 
bulence becomes subsonic and the increase in the internal en- 
ergy due to the shock dissipation ceases, provided the mag- 
netic energy is smaller than the kinetic energy up to that time. 

Therefore, the relativistic turbulence model for GRBs (Lyu- 
tikov 2006; Narayan & Kumar 2009; Kumar & Narayan 2009; 
Lazar et al. 2009) does not seem to work because the rela- 
tivistic turbulence decays much more rapidly than the eddy- 
turnover time fdec ^ feddy ^ / c (sec also Zrake, Jonathan & 
MacFadyen 2010). 

5. SUMMARY OF SIMULATIONS AND COMPARISONS WITH 
PREVIOUS STUDIES 
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Fig. 12. — Top panel: Evolution of the maximum (solid) and average 
(line and points) magnetic field strengths. Middle: Evolution of the magnetic 
(thick solid), kinetic (thick dashed) and internal (thick dotted) energies. Thin 
solid line represents eq (5) with / = 1.8. Botom panel: Evolution of es- 

Using three-dimensional special-relativistic MHD simula- 
tions, we have studied the generation of turbulence as a con- 
sequence of the interactions between a relativistic shock wave 
and density fluctuations. We have found that the magnetic 
field is amplified due to the turbulence dynamo effect and fol- 
lows a power-law decay in later stages. The following items 
are noteworthy: 

• The velocity dispersion of turbulence induced by the 
density fluctuation-shock interaction roughly depends 
linearly on the dispersion of the preshock density inho- 
mogeneity. It can be as large as the postshock sound 
speed when the dispersion of density fluctuations is on 
the order of the mean value (An/n ~ 1). Therefore, 
even if the shock is relativistic, the induced turbulence 
in our situation cannot be as highly relativistic as pos- 
tulated in Narayan & Kumar (2009) model. 

• For transonic and subsonic turbulence, the induced tur- 
bulence maintains its strength for a few eddy-turnover 
times and then decays. In such a turbulent medium, un- 
til the turbulence begins to decay, the magnetic field is 
amplified exponentially in time according to eq. (5) be- 
cause of the effect of field-line stretching. The degree of 



amplification of the magnetic energy before the start of 
the turbulence decay is approximately 100, regardless 
of the initial velocity dispersion of turbulence (eq. [6]). 
The magnetic field continues to grow until the magnetic 
energy becomes comparable to the kinetic energy (eq. 
[7]) after which the magnetic energy follows a power- 
law decay with an exponent of ^ 0.7 (i.e., cx f""^). 

• The evolution of cb (the ratio of the magnetic energy 
density to the internal energy density) in the transonic 
and subsonic turbulence is similar to the magnetic en- 
ergy, because the internal energy is almost constant dur- 
ing the evolution. We have found that when the initial 
eg immediately behind the shock is a few times lO""*, eg 
can grow on the order of 0.1 at maximum (see Fig. 6). 

• The critical length scale below which the back reaction 
of magnetic field on turbulence becomes effective is 
^ 1 /10th the initial inhomogeneity scale. At this scale, 
the magnetic energy becomes comparable to the kinetic 
energy, and the Kolmogorov cascade transitions to the 
MHD critical balance cascade. In addition, below this 
scale, spatial correlations of the local magnetic field ori- 
entation appears (see §4.1.2). 

• For relativistic supersonic turbulence, the kinetic en- 
ergy decay begins an order of magnitude faster than for 
the transonic case, because the formation of a number 
of shock waves directly dissipates the kinetic energy. 
The evolution of eg is completely different from the 
transonic and subsonic cases-it decreases even while 
the magnetic energy is growing because the increase in 
the internal energy by shock dissipations is faster than 
the magnetic field amplification. 

In the remainder of this section, we compare our re- 
sults with previous related studies. Using three-dimensional 
relativistic MHD simulations, Zhang et al. (2009) recently 
showed that the turbulence induced by the Kelvin-Helmholtz 
instability (KHI) in a relativistic shear flow amplifies the mag- 
netic field. They found that eg converges to 5 x 10"^ irre- 
spective of its initial values of 10"'' and 10"^. Our results, 
on the other hand, show that eg depends significantly on the 
initial value. The difference between the two results stems 
from the sources of turbulence. In the simulations of Zhang 
et al. (2009), the shear flow is continuously injected by hand, 
which constantly drives the turbulence via the KHI. On the 
other hand in our simulations, the turbulence that is expected 
to arise from the RMI is given only initially. Considering 
these differences, we find that the results of the two simu- 
lations are essentially the same in the sense that the magnetic 
energy is deposited by turbulent flows induced by a hydro- 
dynamic instability that saturates once it becomes as large as 
kinetic energy, and then the magnetic energy evolves along 
with the kinetic energy (see. Fig. 2 of Zhang et al. 2009 and 
Fig. 5 herein). 

The exponential growth of magnetic energy and the sub- 
sequent power-law decay obtained from our simulations is 
similar to the case of the Weibel instability (Chang et al. 
2008; Keshet et al. 2009). The Weibel instabihty is a pow- 
erful mechanism for the generation and amplification of mag- 
netic fields. However, the typical scale of the magnetic field 
on the order of the plasma skin depth may be insufficient to 
contribute to particle acceleration. In addition, the magnetic 
fields generated at shock front decay rapidly. The advantage 
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of magnetic field amplification by tiie RMI is in its spatial and 
time scales. Because the scale of the RMI is determined by 
that of the preshock density fluctuations, the typical scale of 
magnetic field can be macroscopic and typically comparable 
to the causally connected scale (i.e., the maximum scale in the 
observable region). 

6. DISCUSSION: IMPLICATIONS FOR GAMMA-RAY BURSTS 

The results of our simulations can address various issues. 
For example, we consider here the impUcations for the GRB 
emissions. The initial conditions of our simulations in Table 
1 correspond to the GRB-emission region particularly to the 
internal shocks of GRBs with the following physical quanti- 
ties: 

• The kinetic luminosity is L = Anr^mpnc^T^ ^ 6 x 10^^ 
erg S-' (r/lO'"* cm)2(n/10'° cm'^) (T/lOY, which is a 
typical value, where F is the bulk Lorentz factor. 

• The magnetic energy, carrying a fraction ~ e^^o of 
the total luminosity, is subdonainant. Considering the 
frozen-in magnetic field transported from the central 
engine Bfz, the ratio of the magnetic to the total lumi- 
nosity, Lb,JL = 47rr2(fi|/87r)cr7L = Bl/Sirmpnc^ = 
eso, is nearly conserved during the free expansion of the 
fireball because the area normal to the toroidal magnetic 
field is proportional ^ to oc r~^; hence, the comoving 
toroidal field evolves as Bfz oc r"'r~'. The magnetic 
fraction eB„ corresponds to a central-engine magnetic 
field of 



1/2 



Bo 



/ 87reB„L \ 



'SxlOi^Ge^f 



• We expect a density inhomogeneity in a GRB jet be- 
cause the angular size of a causally connected region 
F"' ~ 10"^ is usually smaller than the jet opening an- 
gle ~ 0.1. The density inhomogeneity is also sug- 
gested by the observations such as the large variation 
in the prompt luminosity compared to that in the af- 
terglow (Kumar & Piran 2000), the spectral and tem- 
poral varieties (loka & Nakamura 2001; Yamazaki et 
al. 2004), and the variabilities of the early afterglow 
(loka et al. 2005) and its polarization (Toma et al. 
2009, and references therein). The comoving size of 
the causally connected region is ^ r/F ^ lO" cm 
(r/lO'"* cm)(F/10-^)"', which may be the typical scale 
of the density fluctuations. The fluctuation scale might 
be smaller than this scale, because the sound velocity is 
less than the light speed c before the shocks. In either 
case, we may take = (^ simulation box size) 
as the fluctuation scale because the simulation is scale 
free. 

In the following, we consider two leading models of GRB 
prompt emission: the synchrotron model and the photosphere 
model. We also apply the afterglow model. 

6.1. Synchrotron Model 

The internal shocks convert kinetic energy into internal en- 
ergy, which goes into the magnetic field and the electron ac- 
celeration with energy fractions cb and e^. The electrons ra- 
diate synchrotron emission that is observed as the prompt 
GRB. This is the internal-shock synchrotron model (Meszaros 
2006). 

The characteristic synchrotron frequency is 



10^^ erg s 



The preshock magnetic fields assumed in our simula- 
tions correspond to an initial fireball at r = ro being 
magnetized as Bo ~ lO^^ G for run Bl, Bq ~ lO" G 
for runs B2, B2-s, and B2-r and B,, 10^" G for run 
B3. In these estimations, we have considered the effect 
of compression by internal shock (e^.ps — Bp^/En ~ 
100 cbq, where the subscript "ps" indicates the values 
of the post internal shock or the initial condition of the 
series of run B). 

• The comoving temperature of the freely expanding fire- 
ball is adiabatically cooled to 

kBT^ksTor~\r/r„Tr^/^ 
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lO^ cm/ 

whereas it rises to ^ 1 GeV after the internal shocks 
(followed by a similar adiabatic cooling). In the simu- 
lation, we initially consider ~ 9.38 MeV to ensure 
the numerical stability. However, this does not affect 
our conclusions as long as the initial temperature is well 
below the postshock temperature. 

' The area normal to the toroidal field evolves oc r'^ (not oc r"') at r > 
roF- ~ 10'^ cm because the shell thickness broadens. However, the shell 
broadening leads to the internal shoclcs between the successive shells, which 
prevent more than twice broadening. 



where fe is the fraction of electrons that are accelerated, TZ is 
the number ratio of electrons (plus positrons) to protons, F is 
the relative Lorentz factor between shells, Af is the variability 
time, 7„ « f~^TZ~^T€eimp/me) is the characteristic random 
Lorentz factor of electrons, and we use 47rr^(B^/87r)cF^ = 
egLint « CBLj/ce and r = 2F^cA/. In the synchrotron model, 
we identify the characteristic frequency with the observed 
peak energy of the Band spectrum (e.g., Zhang & Meszaros 
2002). Our simulations indicate that the RMI can yield a suf- 
ficient magnetic fraction eb ^ 10"^ to reproduce the observed 
peak energy. The necessary condition is that the initial mag- 
netic fraction immediately behind the shock is eB,ps ^ 10 
(see Fig. 6), i.e., the central engine magnetic field is Bq > 10" 
G from Eq. (8). In other words, the Poynting flux may be sub- 
dominant. It is interesting that the maximum level of magnetic 
energy depends on the initial magnetic field in the turbulence 
caused by the RMI instability (see § 4). We emphasize that 
a strong requirement for the amplitude of the initial density 
fluctuations is not necessary because the hundredfold growth 
of magnetic energy can be realized irrespective of the veloc- 
ity dispersion of turbulence (see, eq. [6]). As the dispersion 
of the density fluctuations is reduced, the velocity dispersion 
decreases, which results in a longer timescale for amplifica- 
tion. Because the local magnetic field strength is distributed 
as shown in Fig. 8, the typical magnetic fraction determines 
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the peak energy. Note that the factor f^^TZ'^Teg ^ 10 in 
Eq. (10) may require an efficient electron acceleration ^ 1, 
a large relative Lorentz factor F ~ 10, few positrons 7?. ~ 1, 
and/or a small fraction of accelerated electrons ~ 0.1 (e.g., 
Eichler & Waxman 2005; Toma et al. 2008). 

Electrons should be scattered by disturbed magnetic fields, 
which is required for the first-order Fermi acceleration. The 
relevant wavelength to resonantly scatter electrons with en- 
ergy 7„OTeC^ corresponds to a Larmor radius of 




If this length-scale is attributed to the RMl instability, the ini- 
tial size of the inhomogeneity is required to within a range 

10 to 1.0 times which is much shorter than r/F. How- 
ever, even if the scale of the density inhomogeneities is much 
larger than R,n, the magnetic field fluctuations induced by the 
RMl can scatter the electrons via magnetic mirror reflections. 
Indeed, Beresnyak et al. (2010) recently studied the transport 
of test particles in MHD turbulence, and found effective scat- 
tering of particles by magnetic bottles formed by large-scale 
slow-mode perturbations. 

The internal- shock synchrotron model has several crucial 
problems, one of which is the cooling problem (Meszaros & 
Rees 2000). The cooling time for electrons with the char- 
acteristic Lorentz factor 7^ is usually much shorter than the 
comoving causal time of ^ FAf, so that almost all electrons 
cool down to a Lorentz factor 7^ ~ Gnmec/aTB^TAt. The 
corresponding cooling frequency Vc = Th'y^eB/meC is 




below the characteristic frequency In this case, the low- 
energy spectral index below the peak energy (!/„) becomes 
F,y 3c which contradicts the harder observations a 

RMl turbulence could solve the cooling problem by con- 
tinuously accelerating elections through the stochastic accel- 
eration so-called the second-order Fermi acceleration in non- 
relativistic cases. The quasi-linear theory for election scat- 
tering gives us the scattering timescale as f~{ = (7r/4)/Rf7L, 
where /r is the energy density fraction of the resonant turbu- 
lence to the background magnetic field, and Ql is flie Larmor 
frequency. The simulation results show the non-linear turbu- 
lences i{\B\) ~ A|fi|), and tiie non-resonant scattering may be 
essential as we mentioned before. In spite of those issues, ex- 
trapolating this formula with /r ^ 1, the scattering timescale 
is estimated as ~ fij^^ which is significantly shorter than the 
cooling timescale. Since the mildly relativistic turbulence 
leads to the energy variance per scattering AE /£ ~ Av/c, the 
heating timescale due to the turbulences may be siginificant. 
Asano & Terasawa (2009) demonstrated that the second-order 
Fermi acceleration balances synchrotron cooUng, so that the 
low-energy spectral index becomes as hard as oc u^/^-u^, 
consistent with the observations. In their Monte Carlo simu- 
lations of acceleration, the mean collision time of pitch-angle 
scattering is assumed to be independent of election energy. 
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According to the measurement of the test-particle collision 
frequency in MHD turbulence by Beresnyak et al. (2010), 
the colUsion frequency is independent of energy because the 
particles are mainly scattered by magnetic bottles formed by 
large-scale slow-mode perturbations. Thus we can expect the 
second-order Fermi acceleration in the post-shock turbulent 
region. Asano & Terasawa (2009) also demonstrated that 
if the electron collision frequency decreases with time, the 
accelerated electrons can reproduce the Band spectrum, in- 
cluding its high-energy side. Their assumption is compatible 
with our simulation results, because the turbulence we are dis- 
cussing is time dependent; its magnetic energy decays with 
time at a later stage, which leads to decreasing collision fre- 
quency. The direct numerical simulations including the parti- 
cle acceleration, which would be feasible by applying the test- 
particle approximation, is necessary to quantitatively confirm 
these expectations. 

Another possibility to solve the cooling problem might be 
to consider that the eddy scale (i.e., the density fluctuation 
scale) is much smaller than the comoving causal scale. As 
shown in Fig. 5, the magnetic field decays within a few eddy- 
turnover times. If the decay timescale is comparable to the 
cooling timescale of electrons of 7^, emissions from cooled 
electrons are suppressed. Then, the low-energy specttal index 
becomes that of the synchrotton, oc v^^^, consistent with 
the observations (Pe'er & Zhang 2006). However note that 
the radiative efficiency becomes too small unless the decay 
times is fine tuned to the cooling time. 

We also comment on the jitter radiation (Medvedev 2000), 
in which small-scale turbulences lead to average deflections 
much smaller than the beaming angle and a low-energy spec- 
tral index F^ cx harder than that for the synchrotron, as ob- 
served in a fraction of GRBs. Because the Larmor radius 
Ri ^ jetrieC^ /eB, or more precisely RiHe, is much smaller 
than the typical magnetic field scale, which is roughly one- 
tens of the scale of initial density fluctuations (below which 
the power spectrum decreases as shown in Figs. 7 and 10), 
the jitter radiation does not work for macroscopic RMl turbu- 
lence. 

An interesting prediction of the synchrotron model with 
RMl turbulence is the polarization of the prompt GRB. Since 
the typical scale of the magnetic field is ^ i-/30 (i.e., ^ 
1/lOth the fluctuation scale of ^ the number of do- 

mains with coherent magnetic field would be at least A'^ ^ 10^. 
Therefore, the polarization degree should be less than 

nL<^~2%, (13) 

(Gruzinov & Waxman 1999), which may be probed by the 
future X-ray polarimetric observations (Toma et al. 2009). If 
the polarization of the prompt GRB is above this Umit, we 
have to consider mechanism other than the synchrotron model 
via RMl turbulence, such as the magnetic field advected from 
the central engine. 

Even if the Band component at MeV range is produced by 
a mechanism other than synchrotron (e.g., photosphere emis- 
sion, which is discussed in the next section), the extra high- 
energy component recently identified by the Fermi satellite 
(Abdo et al. 2010, 2009) may require emissions from non- 
thermal particles. The extra components can be explained by 
several models such as early onset of the afterglow (Ghisellini 
et al. 2009; Kumar & Duran 2009), upscattering of exter- 
nal/photospheric photons (Toma et al. 2009b, 2010; Pe'er et 
al. 2010), and hadronic pair cascade (Asano et al. 2009, 2010). 
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For those models, the magnetic field amplification and tur- 
bulences to produce non-thermal particles are indispensable. 
Another interesting method to emit GeV photons is the inter- 
nal shock synchrotron with high Lorentz factor recently pro- 
posed by loka (2010). In such cases the coohng frequency can 
reach fc > 100 GeV if the Lorentz factor is as high as F > 10"* 
in Eq. (12). With eq. (10) {v^ cx F"^), the rising segment of the 
vF^ oc v^~P^I'^ spectrum is stretched down below 1 keV over 
more than 7 energy digits, as observed. In the high Lorentz 
factor model, the maximum synchrotron frequency is lim- 
ited by the magnetic field decay time /^FAf > n'-^enieC / qeB, 
where fs is the ratio of the decay time to the comoving causal 
time and the right hand side is the electron acceleration time 
(with K ~ 1 for the Bohm limit) (loka 2010). This yields 



Bdecay 



100GeV/i«-2(:^) 



104 



\ lO^-' erg s 



3/2 



Af 
10^ 



(14) 



Therefore, to produce high-energy photons, /b ~ 1 is nec- 
essary, i.e., the eddy-turnover time t^y ^ L^l should be 
comparable with the causal time. 

6.2. Photosphere Model 

Another problem with the internal shock synchrotron model 
is an efficiency problem (Kumar 1999; Zhang et al. 2007; 
loka et al. 2006). The observed high-radiative efficiency re- 
quires large dispersion for the Lorentz factor (Kobayashi & 
Sari 2001), which tends to destroy the observed correlations, 

1/2 

L,/ (Yonetoku et al. 2004), because the peak energy v,„ 
is also sensitive to F and F in Eq. (10) (Asano & Kobayashi 
2003). 

The difficulties of the internal shock models lead to the re- 
examination of the original fireball model, in which photons 
are released in a photospheric emission when the fireball be- 
comes optically thin (e.g., Meszaros & Rees 2000; Ramirez- 
Ruiz 2005; Thompson et al. 2007; Ryde & Pe'er 2009; loka 
et al. 2007; loka 2010). The original problem is alleviated 
by introducing dissipation under the photosphere, which can 
bring the thermal peak into the observed range. The photo- 
sphere model can naturally achieve the high efficiency and 
the hard low-energy spectrum. The only crucial flaw is that 
the spectrum tends to be thermal without the nonthermal tails 
observed in GRBs. The electrons with mildly relativistic tem- 
perature can upscatter the thermal photons to produce a high 
energy tail (Beloborodov (2010), see also similar simulations 
for AGNs, Asano & Takahara (2007, 2009)), but flie mecha- 
nism to keep electron temperature higher than photon temper- 
ature is not self-evident. 

A probable solution is the RMl turbulence excited just be- 
low the photosphere, which could convert the thermal spec- 
trum into the observed nonthermal spectrum through Comp- 
tonization. This is similar to the Thompson (1994) model, 
which employs Alfv6n turbulence rather than the RMI tur- 
bulence. If the turbulent velocity is mildly relativistic and 
the Thompson optical depth is about unity (i.e., just below 
the photosphere), the Compton >-parameter is approximately 
unity. In this case, if the turbulent kinetic energy is compa- 
rable or larger than the radiation energy, photons are statisti- 
cally upscattered above the peak energy into a broken power- 
law (Band-like) spectrum. Note that electrons are, in a sense, 
continuously heated because they move together with protons 



in macroscopic turbulence, i.e., « 1 is effectively achieved. 
For the high-radiative efficiency of GRB, it is necessary that 
the turbulent energy is just comparable with the radiation en- 
ergy. In the RMI turbulence model, this is naturally achieved 
by the relativistic shock just below the photosphere, where 
the density is so low that the internal energy (i.e., microscopic 
motion of matter) is not effectively thermalized into radiation 
in the post shock. Hence, the pressure balance between the 
radiation and the matter behind the shock automatically leads 
to near equipartition of radiation and turbulent energy (loka 
2010). 

According to the simulation presented in Fig. 2, the turbu- 
lent velocity dispersion is larger for larger density fluctuations 
audits maximum is Av^ 0.6 x c/\/3 0.3c (run Al). In this 
case, the Compton );-parameter is 



y = TT- 



4m 



1 



Av 
03^ 



(15) 



where the dependence on tt is not square because the fire- 
ball is expanding with decreasing tt- The high-energy spec- 
tral index in cx j/'"^" is given by the unsaturated Compton 
spectrum (e.g., Rybicki & Lightman 2004) with 




-1, 



(16) 



as observed in the density fluctuation An/no ~ 1 (i.e., Av ~ 
0.3c) just below the photosphere tj ^ 5. 

Note that the photosphere model does not explain the high- 
energy photons above ^ TnieC^ ~ 1 GeV (T /Kfi). However 
the high-energy emission can be produced by the subsequent 
emission such as the internal shock and afterglow emission 
(see Sec. 6.1). 

6.3. Afterglow 

The afterglow is thought to be produced by the relativis- 
tic shock between the outflow and the ambient medium via 
synchrotron emission, although the interpretation of early af- 
terglows has not been settled (Zhang et al. 2007; loka et 
al. 2006). The broadband modeUng suggests various mag- 
netic fractions of 10"^ 1^ (b 10~' with a mean of roughly 
cb ~ 10~^ (Panaitescu & Kumar 2002). The magnetic field 
amplification by the compression of ^ /xG circumburst mag- 
netic field merely yields cb ~ 10"', which is too weak for 
the afterglow emission. The small-scale field produced by 
the plasma instabilities such as the Weibel instability decays 
rapidly and does not persist over the emission region (Chang 
et al. 2008; Keshet et al. 2009). 

The RMI turbulence dynamo could be responsible for the 
magnetic field generation of the afterglow. In this mechanism, 
the maximum magnetic energy depends on the initial condi- 
tions (see §4.1.1). If the initial magnetic field is ^ 10"' 
as expected for the compression of a circumburst magnetic 
field, the maximum cb would be <C 10""* because the maxi- 
mum would be smaller than than that for run B3 whose initial 
eg is ^ 10"^. Our maximum magnetic energy level is less 
than that expected in Goodman & Macfadyen (2008), where 
almost all the kinetic energy induced by the density bump- 
shock interaction was supposed to be converted into the mag- 
netic energy. However, other mechanisms such as the cosmic- 
ray or secondary e='= pair streaming instability (Lucek & Bell 
2000; Ramirez-Ruiz et al. 2007) could preamplify the mag- 
netic field to a moderate level es ^10"^, which can be boosted 
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to the necessary level eg ^ 10"^ by the RMI turbulence dy- 
namo. The advantage of this scenario is that preamplification 
may be moderate and not necessarily complete. In addition, 
the dependence of the initial conditions can diversify to the 
magnetic fraction cb as inferred from the observations. 

RMI turbulence could also cause the shallow decay phase, 
which is the most enigmatic feature in the early afterglow. The 
point is that the maximum value of es depends on the initial 
density fluctuation if the the magnetic energy immediately be- 
hind the shock is more than two orders of magnitude smaller 
than the kinetic energy of turbulence (see, eqs. [6] and [7]). In 
addition, depending on An/n and its initial scale Aq, the mag- 
netic field growth timescales fgddy ^ Aq/ Av ^ Aon/ An /cj) 
can be longer than the dynamical timescale R/{cT). This ef- 
fect may also lead to the effective An/n dependence of cb- 
The stellar wind (e.g.. Castor et al. 1975) or the ionizing ra- 
diation (e.g., Bertoldi 1989) from the GRB progenitors may 
diminish the nearby density fluctuations, or induce some hy- 
drodynamical instabilities (Owocki et al. 1988; Ramirez-Ruiz 
et al. 2005). Such radial dependencies of An/n may effec- 
tively cause the evolution of eg (and probably also e^) as dis- 
cussed in Yost et al. (2003) and loka et al. (2006) to explain 
the shallow decay phase of the early afterglow. 

The observed afterglow polarization of 11/, 1% (e.g. 
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